Suppression of elliptic flow in a minimally viscous quark-gluon plasma 
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We compute the time evolution of elliptic flow in non-central relativistic heavy-ion collisions, using 
a (2+l)-dimensional code with longitudinal boost-invariance to simulate viscous fluid dynamics in 
the causal Israel-Stewart formulation. We show that even "minimal" shear viscosity rj/s — ft/(47r) 
leads to a large reduction of elliptic flow compared to ideal fluid dynamics, raising questions about 
the interpretation of recent experimental data from the Relativistic Heavy Ion Collider. 
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PACS numbers: 25.75.-q, 25.75.Ld, 47.75. +f, 12.38.Mh 

The success of the hydrodynamic model in describing 
the bulk of hadron production in Au+Au collisions at 
the Relativistic Heavy Ion Collider (RHIC) Q has led 
to a paradigmatic shift in our view of the quark-gluon 
plasma (QGP): Instead of behaving like a gas of weakly 
interacting quarks and gluons as naively expected 
on the basis of asymptotic freedom in QCD, its collec- 
tive properties rather reflect those of a "perfect liquid" 
with (almost) vanishing viscosity. However, due to quan- 
tum mechanical uncertainty no fluid can have exactly 
zero viscosity Q, and recent work 0] on strongly cou- 
pled gauge field theories, based on techniques exploit- 
ing the AdS/CFT correspondence, suggests an absolute 
lower limit for the ratio of shear viscosity r\ to entropy 
density s: rj/s> h/Air. This raises the interesting ques- 
tion how close to this limit the actual value of the shear 
viscosity of the QGP created at RHIC is. 

Answering this question requires hydrodynamic simu- 
lations for relativistic viscous fluids in which the ratio 
rj/s enters as a parameter. To study the anisotropic ("el- 
liptic") collective flow in non-central heavy-ion collisions, 
from which limits on rj/s can be extracted Q, requires 
a code that evolves the hydrodynamic fields at least in 
the two dimensions transverse to the heavy-ion beam. In 
this Letter we present our first results from such simula- 
tions Q ; a longer paper with a discussion of all technical 
details of our approach is in preparation Q- 

Relativistic hydrodynamics of viscous fluids is techni- 
cally demanding. The straightforward relativistic gen- 
eralization of the non-relativistic Navier-Stokes equation 
yields unstable equations that can lead to acausal sig- 
nal propagation. A causally consistent theoretical frame- 
work was developed 30 years ago by Israel and Stewart 
Q. It involves the simultaneous solution of hydrody- 
namic equations for a generalized energy-momentum ten- 
sor containing viscous pressure contributions, tt^ v (x), to- 
gether with kinetic evolution equations, characterized by 
a (short) microscopic collision time scale t w , for the dy- 
namical approach of towards its Navier-Stokes limit. 
Compared to ideal fluid dynamics, this leads effectively 
to more than a doubling of the number of coupled partial 
differential equations to be solved 



The last couple of years have seen extensive activ- 
ity in implementing the Israel- Stewart equations (and 
slight variations thereof) [1,0, [HI 11, 13] numerically, for 
systems with boost-invariant lon gitu dinal expansion and 
transverse expansion in zero one [ll jl3l FLU Il5j 

and two dimensions 16|, 11 1 (see also Ref. |l8| for a nu- 
merical study of the relativistic Navier-Stokes equation 
in 2+1 dimensions). It is probably fair to say that the 
process of verification and validation of these numerical 
codes is still ongoing: While different initial conditions 
and evolution parameters used by the different groups of 
authors render a direct comparison of their results dif- 
ficult, it seems unlikely that accounting for these differ- 
ences will bring the various published numerical results 
in line with each other. 

We here present results obtained with an indepen- 
dently developed (2+l)-dimensional causal viscous hy- 
drodynamic code, VISH2+1 [i^l^], which has been ex- 
tensively tested (for details see 7|): (i) in the limit of van- 
ishing viscosity, it accurately reproduces results obtained 
with the (2+l)-d ideal fluid code AZHYDRO (ii) for 
homogeneous density distributions (i.e. in the absence of 
density gradients) and vanishing relaxation time it accu- 
rately reproduces the known analytic solution of the rel- 
ativistic Navier-Stokes equation for boost-invariant lon- 
gitudinal expansion (iii) for very short kinetic relax- 
ation times our Israel-Stewart code accurately reproduces 
results from a (2+l)-d relativistic Navier-Stokes code, 
under restrictive conditions where the latter produces nu- 
merically stable solutions; and (iv) for simple analytically 
parametrized anisotropic velocity profiles the numerical 
code correctly computes the velocity shear tensor that 
drives the viscous hydrodynamic effects. 

VISH2+1 solves the equations for local energy- 
momentum conservation, d m T mn = 0, with 



T'' 



eu m u n -pA mn + ir mn , A mn = g mn -u m u n , (1) 



together with kinetic equations for the viscous shear pres- 



sure 7T 



D7t mn = , 2 mn^mm _ ( u m Tr nk +u n Tr mk )Du k . (2) 
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Here D = u m d m is the time derivative in the local co- 
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moving frame and a r 



±A mn d, 
and traceless 



V (m u n> = ±(V"V' l +V n w m ) - 
(with W m — A ml di) is the symmetric 
velocity shear tensor. We use a 
fixed Eulerian space-time grid in curvilinear coordi- 
nates x m = (t, x, y, 77), with longitudinal proper time 
t = \/t 2 —z 2 and space-time rapidity r\ = | In |i£ where 
z is the beam direction and (x, y) are the two trans- 
verse directions. d m indicates the covariant deriva- 
tive in direction x m in this coordinate system. As in 
Refs. " 



u,\m\m, u 



|18| | we neglect bulk viscosity and 
heat conduction as presumably subdominant effects in a 
QGP with approximately vanishing net baryon density. 

We implement longitudinal boost-invariance via the 
ansatz u m = (u T ,u x ,u v , u v ) = j±(l, v x , v y ,0), with j± = 



-1/2 



using ^-independent initial conditions. 



The equations to be solved are 3 hydrodynamic equations 
for T TT , T TX and T Ty , together with 4 kinetic equations 
for tt^, 7r TT , n TX and w TV 0, E2]. To check the numer- 
ics we also evolved additional, redundant components of 
TT mn and confirmed that the identities u m ir mn = = 7r,™ 
are preserved over time. 

Due to the limited size of the transverse {x, y) grid in 
our current code, we here present results only for Cu+Cu 
collisions; simulations of the larger Au+Au system will 
soon be forthcoming. We use standard Glauber model 
initial conditions [H, assuming wounded- nucleon scaling 
of the initial transverse energy density e(x,y,To;b) (23[, 
with a Woods-Saxon radius Rc u — 4.2 fm, surface thick- 
ness £ = 0.596 fm, and equilibrium density po = 0.17fm -3 . 
We scale this profile to a peak initial energy density 
in central {b = 0) collisions of eo = 30 GcV/fm 3 . This 
is higher than expected for Cu+Cu collisions at RHIC 
but ensures that the system spends enough time in the 
QGP phase to explore the effects of shear viscosity on 
the evolution of anisotropic flow in this phase. We start 
the hydrodynamic evolution at To = 0.6fm/c with van- 
ishing transverse flow, both for the viscous evolution and 
ideal fluid dynamical comparison runs. For the equa- 
tion of state (EoS) we use a slight variation of EOS Q 
from Ref. [23| (which implements a phase transition at 
T c = 164MeV between a free quark-gluon gas above T c 
and a chemically equilibrated hadron resonance gas be- 
low T c ) where the sharp corners at either end of the 
Maxwell contraction at T c are rounded off for numeri- 



cal stability [24(. Neither the initial conditions nor the 
EoS have been fine-tuned for a realistic comparison with 
experimental data; we here emphasize the comparison 
between ideal fluid and viscous hydrodynamic evolution, 
in order to identify qualitative differences between the 
two and to quantitatively understand their origin. 

The viscous hydrodynamic equations contain two pa- 
rameters, the shear viscosity 77 and the kinetic relaxation 
time tv. All simulations presented here are done with 
"minimal" viscosity Q rj/s — 1/4tt, while tv is varied 
between \ (default) and \ of the classical Boltzmann gas 
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FIG. 1: (Color online) Time evolution of the central entropy 
density (a) and average radial flow velocity (b) in central 
Cu+Cu collisions (see text for details). Here and later stars 
indicate the time when all matter is frozen out. 



estimate rf tz = ^ = ^^U>M- 

In Figure [T] we show the evolution of the central en- 
tropy density s(r=0) and the average radial flow (v r ) 
(with the Lorentz contracted energy density j±e as 
weight function) for central (b=0) Cu+Cu collisions. The 
curves labeled "(0+l)-d" correspond to I-dimensional 
boost-invariant longitudinal expansion without trans- 
verse flow (i.e. to transversally homogeneous initial con- 
ditions). One sees that in this case (which, for a simple 
EoS and in the Navier-Stokes limit, can be solved ana- 
lytically 0]) shear viscosity reduces the cooling rate, due 
to a well-known reduction of the work done by longitudi- 
nal pressure Q - additional entropy is produced, and the 
entropy density decreases more slowly than the I/r-law 
for ideal fluids. For transversally inhomogeneous initial 
conditions ("(l+I)-d hydro"), the developing radial flow 
increases the cooling rate compared to the case without 
transverse expansion, for both ideal and viscous fluid dy- 
namics. However, in this case the shear viscous effects 
increase the transverse pressure relative to the ideal fluid 
case, leading to a faster build-up of radial flow (bottom 
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panel in Fig. [I]), which in turn accelerates the cooling of 
the fireball center in such a way that by the time the ex- 
pansion becomes effectively three-dimensional (indicated 
by the r -3 line in Fig. QIa)) the viscous fluid cools more 
rapidly than the ideal one (This was also seen (al- 
though not emphasized) in Rcfs. [n|, Q EEL The 
larger radial flow of the viscous fluid (Fig [lib)) leads to 
flatter final transverse momentum spectra [14U15L 
but (as pointed out in [H, this effect can be largely 
compensated by starting the viscous hydrodynamic evo- 
lution later and with lower initial energy density. 
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FIG. 2: (Color online) Time evolution of the dominant com- 
ponents of the shear viscous pressure tensor, normalized by 
e+p and averaged over the transverse plane, for two different 
initial conditions (red and green, respectively). Note that the 
normalization factor e+p ~ T 4 decreases rapidly with time. 



The dotted and dash-dotted lines in Fig.QTb) show that 
(at least for the short relaxation time r T = ir^ oltz = 
0.24 (200McV) fm / c cons idered here) the initial condi- 
tions for the shear viscous pressure tensor don't matter 
much: whether 7r m ™ is initially taken to vanish or to as- 
sume its Navier-Stokes limit 2r/a mn , it quickly relaxes 
to the same function at times t—tq > 4tv « lfm/c. 
This is seen more explicitly in Fig [2] where we show, for 
non-central Cu+Cu collisions at 6 = 7fm and the same 
two sets of initial conditions for the viscous shear pres- 
sure tensor, the time evolution of the dominant compo- 
nents of Tr mn , normalized to the equilibrium enthalpy 
e+p (which sets the scale for the energy-momentum ten- 
sor in the ideal fluid limit) and averaged over the trans- 
verse plane. Other components of (ir mn ) are at least 
an order of magnitude smaller than the ones shown. 
The signs of (tt"") < and (Tr xx ),(w yy ) > reflect 
the reduced longitudinal and increased transverse pres- 
sure caused by shear viscosity. The negative difference 
((tt xx — n yv ) / (e+p)) < seen in Fig.[2]causes a significant 
viscous reduction of the total momentum anisotropy e p 
which we discuss next. 

The anisotropic flow is driven by the spatial source 
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FIG. 3: (Color online) Time evolution of the spatial eccen- 
tricity e x (a) and momentum anisotropy e p (b) for non-central 
Cu+Cu collisions at b— 7fm. See text for details. 



eccentricity e x 



(where 



denotes the en- 



_ (y 2 +x 2 ) 

ergy density weighted average over the transverse plane 
[lj]) and the anisotropic pressure gradients it generates. 
Fig. [3K shows that it decreases as a function of time, and 
does so more rapidly initially, but more slowly later for 
the viscous fluid than in ideal hydrodynamics. The ini- 
tial drop rate for e x depends on the initial value for the 
viscous pressure tensor, but after about 1 fm/c different 
initializations for 7r mn lead to parallel evolution histo- 
ries for e x . The largest initial decrease for e x is observed 
for the largest initial viscous pressure tensor; initial free- 
streaming of the matter would correspond to an extreme 
case of viscous fluid dynamics, leading to even larger ini- 
tial TT mn and even faster initial decrease of the spatial 
eccentricity than shown in Fig. [3K [Isj]- 

In Fig. [3j3 we plot, for ideal and viscous hydrody- 
namic evolution, the total momentum anisotropy av- 
eraged over the transverse plane - — — 

(T XX +7r IX —Tq V —7r yy ) 



I +Tg y +nyy) ' 



P (r** +Tyy) 
as a function of time. The closely 



spaced green and red lines distinguish different initial- 
izations for TT mn , as specified in Fig. [3^, showing weak 
sensitivity to these initial values. More interesting is the 
separation of the contributions to e p arising from the ideal 
fluid part T™" (which only tracks the differences in the 
evolution of flow velocity and thermal pressure between 
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ideal and viscous hydrodynamics) and from the viscous 
pressure components 7r m ™ (dashed and dotted lines in 
Fig. [3Jd). Initially, the viscous flow and thermal pres- 
sure evolution is not very different from the ideal fluid 
case, and it takes a while until the viscous pressure effects 
manifest themselves in a significant reduction of the flow 
anisotropy, thereby modifying the ideal fluid part J 1 ™" of 
the energy momentum tensor. However, the viscous pres- 
sure components 7r m ™ contribute themselves a negative, 
initially large part to the total momentum anisotropy, re- 
sulting (for 2 — _L) in an overall reduction of the latter 
by almost 50% relative to the ideal fluid case over the 
entire time history. At late times, the viscous pressure 
components become small, but by then their negative ef- 
fect on the buildup of elliptic flow has fully manifested 
itself in the collective flow profile and is thus carried by 
the ideal fluid part J"™" of the energy momentum tensor. 
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FIG. 4: (Color online) (a) The elliptic flow v^ijpr) for pi- 
ons, kaons and protons from ideal fluid dynamics (solid lines) 
and viscous hydrodynamics (dotted and dashed lines). Dot- 
ted lines account only for viscous effects on the flow pattern 
that enters the equilibrium part of the distribution function; 
dashed lines additionally include viscous (non-equilibrium) 
corrections to the latter, (b) Effects of different choices for 
the kinetic relaxation time tv and for the initial viscous pres- 
sure 7r mn on pion elliptic flow (with the same separation of 
equilibrium and non-equilibrium contributions as in part (a)). 



The elliptic flow v£ — (cos(2^ p ))j, calculated as the 



distribution -, , , — . 

dyp T dp T d<fi p ■ 

and measures the distribution of the total momentum 
anisotropy e p over the various hadron species and over 
transverse momentum px- This distribution depends on 
the chemical composition of the fireball at kinetic freeze- 
out, here assumed (unrealistically) to be a chemically 
equilibrated hadron resonance gas. The hadron momen- 
tum distributions are calculated as a Cooper-Frye inte- 
gral [H over a decoupling surface S of constant temper- 
ature T dec = 130MeV: 



dNi 



dyprdprdc/jp 
p ■ d 3 a(x) 
(2tt) 3 



p ■ d 3 a(x) 

: (2^) 3 



f$(x,p)+6f^(x,p) 



(3) 



1 + (l-f {i) (x y))- pmpn ^"fo) 



Here /iq y P t(x) ) IS * ne l° ca l thermal equilibrium dis- 
tribution function with temperature T(x), boosted to 
the laboratory frame by the local hydrodynamic flow ve- 
locity u m (x) - both taken from the hydrodynamic out- 
put on the freeze-out surface £ (whose normal vector is 
d 3 a(x)). Sf^ describes the deviation from local ther- 
mal equilibrium due to viscous effects and is given by 
the last term in the bottom line of Eq. ([3]) @, [l2j]. It 
is proportional to the viscous pressure tensor 7r mn , and 
even though (in contrast to early times) TT mn (x) is small 
on the freeze-out surface, its effect on the distribution 
function grows quadratically with momentum, leading to 
a breakdown of the (viscous) hydrodynamic calculation 
of hadron spectra at sufficiently large pt (in our case 
\6N w (p)/N£ q (p) \ > 50% for p T > 2.5 GeV/c where N*(p) 
denotes the pion momentum spectrum). 

In Fig. [4] we compare, for a few common hadron 
species, the elliptic flow V2{pt) from ideal and viscous 
fluid dynamics. One sees that even "minimal" shear 
viscosity ^ = j- causes a dramatic suppression of ellip- 
tic flow. The effects seen in Fig. [4] seem to be even 
larger than those reported in [l7( - this discrepancy 
clearly calls for clarification [30j . Even without account- 
ing for the slight deviations Sf of the distribution func- 
tion from its thermal equilibrium form, caused by small 
but non-vanishing shear pressure tensor components on 
the freeze-out surface, we see that elliptic flow is reduced 
by almost 50% (dotted lines in Fig. [4|), due to the al- 
ready mentioned reduction of azimuthal anisotropies of 
the hydrodynamic flow field. This effect was not even 
considered by Teaney in Q when he attempted to con- 
strain the shear viscosity of the QGP using RHIC data. 
Teaney's argument was entirely based on the viscous cor- 
rections Sf ~ p m p n "n mn ~ Pt arising from non-vanishing 
viscous pressure on the freeze-out surface (the difference 
between the dotted and dashed lines in Fig. 0]) which, at 
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low pt, is a much smaller effect. His phenomenological 
limit [5| on the shear viscosity of the fireball matter cre- 
ated at RHIC is therefore not restrictive enough: even a 
superficial comparison of the shapes of the i>2(pt) curves 
in Fig. 2] with experimental data [l| suggests that RHIC 
data may be inconsistent with the conjectured lower limit 
2 = j- for the QGP shear viscosity. 

This conclusion, even though tentative since it is not 
yet based on a quantitative data comparison with calcu- 
lations that use a more realistic EoS and better initial 
conditions, appears to be robust since neither a 100% 
variation of the initial value for 7r mn nor a 50% reduction 
of the kinetic relaxation time for the viscous pressure ten- 
sor (see Fig. 0b) are able to significantly attenuate the 
strong viscous reduction of elliptic flow that we see in 
our calculations. While it supports the new paradigm 
of the "perfect fluidity" of the QGP created at RHIC, a 
possible violation of the conjectured "minimum viscosity 
bound" Q is a serious matter that raises the question 
whether other explanations might be possible. Among 
the possibilities that one might contemplate are that the 
initial spatial source eccentricity in RHIC collisions has 
been significantly underestimated [29( or that a realistic 
equation of state is effectively stiffer than the one used 
here. T. Cohen and collaborators 31] have advanced 
counter examples of theories that appear to contradict 
the existence of a universal "minimum viscosity bound" , 
although these examples to not include QCD. Our find- 
ings suggest that QCD might belong to the list of ex- 
ceptions. Lublinsky and Shuryak [32| have shown that 
higher order corrections to the Israel-Stewart theory of 
viscous relativistic hydrodynamics tend to decrease vis- 
cous entropy production, so there may be a possibility 
that they also reduce viscous effects on the elliptic flow 
[30j . If this turns out to be case, our tentative conclusion 
as formulated above is premature. Obviously all these 
issues must be carefully investigated before a complete 
understanding of the RHIC data can be achieved. 
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